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Abstract — We  present  the  Bayesian  consensus  filter  (BCF)  for 
tracking  a  moving  target  using  a  networked  group  of  sensing 
agents  and  achieving  consensus  on  the  best  estimate  of  the  prob¬ 
ability  distributions  of  the  target’s  states.  Our  BCF  framework 
can  incorporate  nonlinear  target  dynamic  models,  heterogeneous 
nonlinear  measurement  models,  non-Gaussian  uncertainties,  and 
higher-order  moments  of  the  locally  estimated  posterior  proba¬ 
bility  distribution  of  the  target’s  states  obtained  using  Bayesian 
filters.  If  the  agents  combine  their  estimated  posterior  probability 
distributions  using  a  logarithmic  opinion  pool,  then  the  sum  of 
Kullback-Leibler  divergences  between  the  consensual  probability 
distribution  and  the  local  posterior  probability  distributions  is 
minimized.  Rigorous  stability  and  convergence  results  for  the 
proposed  BCF  algorithm  with  single  or  multiple  consensus  loops 
are  presented.  Communication  of  probability  distributions  and 
computational  methods  for  implementing  the  BCF  algorithm  are 
discussed  along  with  a  numerical  example. 

I.  Introduction 

Distributed  estimation,  where  a  group  of  networked  agents 
collectively  estimate  the  target’s  states,  can  be  used  for  envi¬ 
ronment  and  pollution  monitoring,  tracking  dust  or  volcanic 
ash  clouds,  tracking  orbital  debris  or  asteroids  in  space,  etc 

[I] -  [5].  The  term  distributed  estimation  refers  to  finding  the 
best  estimate  of  the  target’s  states  using  the  sensor  network 
while  the  term  consensus  means  reaching  an  agreement  across 
the  network  [6]-  [10]. 

Many  existing  algorithms  for  distributed  estimation  [1]-  [5], 

[II] -  [15]  aim  to  reach  an  agreement  across  the  network  on 
the  estimated  mean  (first  moment  of  the  estimated  probability 
distribution)  of  the  target  dynamics,  but  cannot  incorporate 
nonlinear  target  dynamics,  heterogeneous  nonlinear  measure¬ 
ment  models,  non-Gaussian  uncertainties,  or  higher-order  mo¬ 
ments  of  the  locally  estimated  posterior  probability  distribution 
of  the  target’s  states.  It  is  difficult  to  recursively  combine 
local  mean  and  covariance  estimates  using  a  linear  consensus 
algorithm  because  the  dimension  of  the  vector  transmitted 
by  each  agent  increases  linearly  with  time  due  to  correlated 
process  noise  [16]  and  the  covariance  update  equation  is 
usually  approximated  by  a  consensus  gain  [17], 

Multi-agent  tracking  or  sensing  networks  are  deployed  in  a 
distributed  fashion  when  the  target  dynamics  have  complex 
temporal  and  spatial  variations.  Hence,  it  is  necessary  to 
preserve  the  complete  information  captured  in  the  locally 
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estimated  posterior  probability  distribution  of  the  target’s  states 
while  achieving  consensus  across  the  network.  The  main  aim 
of  this  paper  is  to  extend  the  scope  of  distributed  estimation 
algorithms  to  track  targets  with  general  nonlinear  dynamic 
models  with  stochastic  uncertainties,  thereby  addressing  the 
aforementioned  shortcomings. 

Bayesian  filters  [18],  [19]  recursively  calculate  the  probabil¬ 
ity  density/mass  function  of  the  beliefs  and  update  them  based 
on  new  measurements.  The  main  advantage  of  Bayesian  filters 
over  Kalman  filter-based  methods  for  estimation  of  nonlinear 
target  dynamic  models  is  that  no  approximation  is  needed 
during  the  filtering  process.  Advances  in  computational  capa¬ 
bility  have  facilitated  the  implementation  of  Bayesian  filters 
for  robotic  localization  and  mapping  [20]  as  well  as  planning 
and  control  [21].  Practical  implementation  of  these  algorithms, 
in  their  most  general  form,  is  achieved  using  particle  filtering 
[22]  and  Bayesian  programming  [23],  This  paper  focuses  on 
developing  a  consensus  framework  for  distributed  Bayesian 
filters. 

The  statistics  literature  deals  with  the  problem  of  reaching  a 
consensus  among  individuals  in  a  complete  graph,  where  each 
individual’s  opinion  is  represented  as  a  probability  distribution 
[24],  [25];  and  under  select  conditions,  it  is  shown  that  con¬ 
sensus  is  achieved  within  the  group  [26].  Exchange  of  beliefs 
in  decentralized  systems,  under  communication  constraints,  is 
considered  in  [27],  [28].  Algorithms  for  combining  probability 
distributions  within  the  exponential  family,  i.e.,  the  limited 
class  of  unimodal  distributions  that  can  be  expressed  as  an 
exponential  function,  are  studied  in  [29],  [30].  If  the  target’s 
states  are  discrete  random  variables,  then  the  local  estimates 
can  be  combined  using  a  tree-search  algorithm  [31]  or  a  linear 
consensus  algorithm  [32],  [33],  In  contrast,  this  paper  focuses 
on  developing  generalized  Bayesian  consensus  algorithms  with 
rigorous  convergence  analysis  for  achieving  consensus  across 
the  network  without  any  assumption  on  the  shape  of  local  prior 
or  posterior  probability  distributions.  The  proposed  distributed 
estimation  using  Bayesian  consensus  filtering  aims  to  reach 
an  agreement  across  the  network  on  the  best  estimate,  in 
information  theoretic  sense,  of  the  probability  distribution  of 
the  target’s  states. 

In  this  paper,  we  assume  that  agents  generate  their  local 
estimate  of  the  posterior  probability  distribution  of  the  tar¬ 
get’s  states  using  Bayesian  filters  with/without  measurement 
exchange  with  neighbors.  Then,  we  develop  algorithms  for 
combining  these  local  estimates,  using  the  logarithmic  opinion 
pool  (LogOP),  to  generate  the  consensual  estimate  of  the  prob¬ 
ability  distribution  of  the  target’s  states  across  the  network. 
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Finally,  we  introduce  the  Bayesian  consensus  filter  (BCF), 
where  the  local  prior  estimates  of  the  target’s  states  are  first 
updated  and  the  local  posterior  probability  distributions  are 
recursively  combined  during  the  consensus  stage,  so  that  the 
agents  can  estimate  the  consensual  probability  distribution  of 
the  target’s  states  while  simultaneously  maintaining  consensus 
across  the  network.  The  pseudo-code  for  the  algorithm  is  given 
in  Algorithm  1.  The  novel  features  of  the  BCF  algorithm  are: 

•  The  algorithm  can  be  used  to  track  targets  with  general 
nonlinear  time-varying  target  dynamic  models. 

•  The  algorithm  can  be  used  by  a  SC  balanced  network 
of  heterogeneous  agents,  with  general  nonlinear  time- 
varying  measurement  models. 

•  The  algorithm  achieves  global  exponential  convergence 
across  the  network  to  the  consensual  probability  distribu¬ 
tion  of  the  target’s  states. 

•  The  consensual  probability  distribution  is  the  best  es¬ 
timate,  in  the  information  theoretic  sense  because  it 
minimizes  the  sum  of  KL  divergences  with  the  locally 
estimated  posterior  probability  distributions.  If  a  central 
agent  receives  all  the  local  posterior  probability  distri¬ 
butions  and  is  tasked  to  find  the  best  estimate  in  the 
information  theoretic  sense,  then  it  would  also  yield 
the  same  consensual  probability  distribution.  Hence,  we 
claim  to  have  achieved  distributed  estimation  using  the 
BCF  algorithm. 

The  Hierarchical  BCF  algorithms  is  used  when  some  of  the 
agents  do  not  observe  the  target.  We  apply  the  Hierarchical 
BCF  algorithm  to  the  problem  of  tracking  orbital  debris  in 
space  using  the  space  surveillance  network  on  Earth. 

A.  Notation 

The  time  index  is  denoted  by  a  right  subscript  while  the 
agent  index  is  denoted  by  a  lower-case  right  superscript.  Let 
Xk  represent  the  true  target’s  states  at  the  fcth  time  instant.  Let 
zi  represents  the  measurement  taken  by  the  jth  agent  at  the 
/,:lh  time  instant.  Let  Tk  represents  the  estimated  probability 
density  function  (pdf)  of  the  target’s  states  over  the  state  space 
X.  The  symbol  p(-)  also  refers  to  pdf  over  the  state  space  X. 

The  graph  represents  the  directed  time-varying  commu¬ 
nication  network  topology  at  the  kth  time  instant,  where  all 
the  agents  form  the  set  of  vertices  V  and  £k  is  the  set  of 
directed  edges.  Let  A fk  denote  the  neighbors  of  the  fh  agent 
at  the  kth  time  instant  from  which  it  receives  information  and 
J!  :=  Ml  U  {j}  denote  the  set  of  inclusive  neighbors. 

Let  N,  R,  and  Rmx”  be  the  sets  of  natural  numbers  (positive 
integers),  real  numbers,  and  m  by  n  matrices.  Let  A  and  er 
represent  the  eigenvalue  and  the  singular  value  of  a  square 
matrix.  Let  1  =  [1, 1, . . . ,  1]T,  I,  and  0  be  the  ones  vector, 
the  identity  matrix,  and  the  zero  matrix  of  appropriate  sizes. 
The  symbols  |-|,  [•],  ln(-),  and  logc(-)  represent  the  absolute 
value,  ceiling  function,  the  natural  logarithm  and  the  logarithm 
to  the  base  c.  Finally,  ||  •  ||^  represents  the  tp  vector  norm.  The 
£p  function  denotes  the  set  of  all  functions  f(x):  R”x  —>  R 
with  the  bounded  integral  (fx  \f(x)\pd/j,(x))  1^P,  where  /r  is 
a  measure  on  SC. 


II.  Preliminaries 

In  this  section,  we  first  state  four  assumptions  used  through¬ 
out  this  paper  and  then  introduce  the  problem  statement  of 
BCF.  Next,  we  discuss  an  extension  of  the  Bayesian  filter  to 
sensor  fusion  over  a  network. 

Assumption  1.  In  this  paper,  all  the  algorithms  are  presented 
in  discrete  time.  □ 

Assumption  2.  The  state  space  (X  C  R":c)  is  closed  and 
bounded,  hence  X  is  compact.  □ 

Assumption  3.  All  continuous  probability  distributions  are 
upper-bounded  by  some  large  value  Ml  €  R.  □ 

Assumption  4.  The  inter-agent  communication  time  scale  is 
much  faster  than  the  tracking/estimation  time  scale.  □ 

These  assumptions  have  been  introduced  to  make  the  algo¬ 
rithms  computationally  tractable  and  to  take  advantage  of  the 
results  dealing  with  bounded  functions  and  compact  support. 
Note  that  discrete,  continuous  or  mixed  probability  distribution 
can  be  handled  in  a  unified  manner  using  measures.  Hence, 
every  probability  distribution  in  this  paper  is  expressed  as  a 
probability  density  function  (pdf)  over  X. 

Let  X  C  R"x  be  the  n^-dimensional  state  space  of  the 
target.  The  dynamics  of  the  target  in  discrete  time  {xk,k  £ 
N,  Xk  £  Xj  is  given  by: 

xk  =  fk{xk- l,Vk_i),  (1) 

where  fk  :  Rnx  x  R"x  -A  R”x  is  a  possibly  nonlinear  time- 
varying  function  of  the  state  Xk-i  and  an  independent  and 
identically  distributed  (i.i.d.)  process  noise  Vk-i,  where  nv  is 
the  dimension  of  the  process  noise  vector.  Let  m  heteroge¬ 
neous  agents  simultaneously  track  this  target  and  estimate  the 
pdf  of  the  target’s  states  (where  m  does  not  change  with  time). 
The  measurement  model  of  the  jth  agent  is  given  by: 

zl  =  hUxk>wk)>  Vj  G  (2) 

where  hj,  :  R”x  x  Rn”j  — >•  Rra*j  is  a  possibly  nonlinear  time- 
varying  function  of  the  state  Xi-  and  an  i.i.d.  measurement 
noise  wi,  where  nZj,nwj  are  dimensions  of  the  measurement 
and  measurement  noise  vectors  respectively.  Note  that  the 
measurement  model  of  agents  is  quite  general  since  it  ac¬ 
commodates  heterogeneous  sensors,  partial  state  observation, 
varied  sensors  for  sensor  fusion  in  frequency  domain,  etc. 

The  objective  of  the  BCF  is  to  estimate  the  target’s  states 
and  maintain  consensus  across  the  network.  This  objective  is 
achieved  in  two  steps:  (i)  each  agent  locally  estimates  the  pdf 
of  the  target’s  states  using  a  Bayesian  filter,  and  (ii)  each 
agent’s  local  estimate  converges  to  a  global  estimate  during 
the  consensus  stage. 

A.  Bayesian  Filter  with  Measurement  Exchange 

The  objective  of  Bayesian  filtering  with/without  measure¬ 
ment  exchange  is  to  estimate  the  posterior  pdf  of  the  tar¬ 
get’s  states  at  the  kth  time  instant,  which  is  denoted  by 
Ml,\/j  £  { 1, . . . ,  to},  using  the  estimated  prior  pdf  of  the 
target’s  states  Tl_x  from  the  (fc—  l)th  time  instant  and  the  new 
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measurement  array  obtained  at  the  /,:lh  time  instant.  Exchange 
of  measurements  is  optional  since  heterogeneous  agents,  with 
different  priors,  fields  of  view,  resolutions,  tolerances,  etc., 
may  not  be  able  to  combine  measurements  from  other  agents. 
For  example,  if  a  satellite  in  space  and  a  low  flying  quadrotor 
are  observing  the  same  target,  then  they  cannot  exchange 
measurements  due  to  their  different  fields  of  view. 

If  an  agent  can  combine  measurements  from  another  neigh¬ 
boring  agent  during  its  update  stage,  then  we  call  them  mea¬ 
surement  neighbors.  In  this  section,  we  extend  the  Bayesian 
filter  by  assuming  that  each  agent  transmits  its  measurements 
to  other  agents  in  the  network,  and  receives  the  measurements 

S3 

from  its  measurement  neighbors.  Let  zkk  :=  {zk,  W  €  Sk} 
denote  the  array  of  measurements  taken  at  the  kth  time 
instant  by  the  measurement  neighbors  of  the  jth  agent,  where 
Si  C  Ji  denotes  the  set  of  measurement  neighbors  among 
the  inclusive  neighbors  of  the  jth  agent. 

It  is  assumed  that  each  agent  either  knows  the  prior  Fq  = 
p(x0\zJ0)  or  J-Q  is  assumed  to  be  uniformly  distributed  over 
X.  The  prediction  stage  involves  using  the  target  dynamics 
model  (1)  to  obtain  the  estimated  pdf  of  the  target’s  states  at 
the  klh  time  instant  via  the  Chapman-Kolmogorov  equation: 

Pi(xk)  =  Pk(xk\xk-i) Pi-^Xk-i)  dp(x)  (3) 
J  x 


using  Bayesian  filters  with/without  measurement  exchange. 
During  each  of  the  nioop  iterations  within  the  consensus  stage 
in  Algorithm  1,  this  estimate  is  updated  as  follows: 

rt,v=T  (ytej£{rt,v-i})  €  {1,  •  •  • , in], Vis  £  N,  (5) 

where  T(-)  is  the  linear  or  logarithmic  opinion  pool  for 
combining  the  pdf  estimates.  Note  that  the  problem  of  mea¬ 
surement  neighbors  does  not  arise  here  since  all  pdfs  are 
expressed  over  the  complete  state  space  X. 

Let  T\ . . .  . ,  limrwoc  Fn,  F*  be  real-valued  measurable 
functions  on  X,  SC  be  the  Borel  cr-algebra  of  X,  and  sV  be 
any  set  in  SC .  Let  /i  jrn ,  //jr»  denote  the  respective  induced 
measures  of  Fn ,  F*  on  SC . 

Lemma  1.  (Pointwise  convergence  implies  convergence  in 
Total  Variation)  If  Fn  converges  to  IF*  pointwise,  i.e., 
limn^oo  Fn  =  F*  pointwise;  then  the  measure  pj:n  converges 

T.V. 

in  TV  to  the  measure  pj^*,  i.e.,  lim-n^oo  pj —A  pr*. 

Proof:  The  proof  follows  from  Scheffe’s  theorem  and 
dominated  convergence  theorem  [36,  pp.  84].  ■ 

The  first  method  of  combining  the  estimates  is  motivated  by 
the  linear  consensus  algorithms  widely  studied  in  the  literature 
[6]-  [10].  The  pdfs  are  combined  using  the  Linear  Opinion 
Pool  (LinOP)  of  probability  measures  [24],  [25]: 


The  probabilistic  model  of  the  state  evolution  pJk(xk\xk-i) 
is  defined  by  the  target  dynamics  model  (1)  and  the  known 
statistics  of  the  i.i.d.  process  noise  vk-i-  The  new  measure- 

s j 

ment  array  (zF)  is  used  to  compute  the  posterior  pdf  of  the 

gj 

target’s  states  ( Fk  =  p,k(xk\zkk))  during  the  update  stage 
using  Bayes’  rule  (4): 


Pk(xk  IV)  = 


(UeeslPk(zk\xk))  P{(xk) 

lx  ( nfe5£  Pfc(**i*fc))pfc(*fc)  Mx) 


(4) 


H,v  =  E  W-V  V7  e  {!,  •  •  •  £  N,  (6) 

where  a?ku_l  =  1  and  the  updated  pdf  F°kv  after 

the  0th  consensus  loop  is  a  weighted  average  of  the  pdfs  of 
the  inclusive  neighbors  Fk  l/_1,  VI  £  Jk  from  the  (y  —  l)th 
consensus  loop,  at  the  fcth  time  instant.  The  LinOP  solution  is 
typically  multimodal  and  depends  on  the  assumption  that  the 
same  0-1  scale  is  used  by  every  agent,  so  no  clear  choice  for 
jointly  preferred  estimate  emerges  from  it  [25], 


The  likelihood  function  pk{zk \xk),  W  £  Sk  is  defined  by 
the  measurement  model  (2),  and  the  corresponding  known 
statistics  of  the  i.i.d.  measurement  noise  wk. 

Note  that  (4)  is  similar  to  the  empirical  equation  for 
Independent  Likelihood  Pool  given  in  [34]  and  a  generalization 
of  the  Distributed  Sequential  Bayesian  Estimation  Algorithm 
given  in  [35].  The  structure  of  (4)  ensures  that  an  arbitrary  part 
of  the  prior  distribution  does  not  dominate  the  measurements. 

III.  Combining  Probability  Distributions 

In  this  section,  we  present  algorithms  for  achieving  con¬ 
sensus  in  probability  distributions  across  the  network.  The 
objective  of  the  consensus  stage  in  Algorithm  1  is  to  guarantee 
pointwise  convergence  of  each  Fl  to  a  consensual  pdf  Ff 
which  is  independent  of  j.  This  is  achieved  by  each  agent 
recursively  transmitting  its  estimated  pdf  of  the  target’s  states 
to  other  agents,  receiving  estimates  of  its  neighboring  agents, 
and  updating  its  estimate  of  the  target.  Let  FJk  0  =  F3k 
represent  the  local  estimated  posterior  pdf  of  the  target’s  states, 
by  the  jlh  agent  at  the  start  of  the  consensus  stage,  obtained 


A.  Consensus  using  the  Logarithmic  Opinion  Pool 

Note  that  Fk  v  =  jPk  v{xk),fxk  £  X  represents  the  pdf  of 
the  estimated  target’s  states  by  the  jth  agent  during  the  uth 
consensus  loop  at  the  /,:th  time  instant.  The  LogOP  is  given  as 
[37]: 


=P’kAxk ) 


Vj  G  {1, . . . ,  m},Vu  £  N, 


(7) 


where  Yleejj  aku-i  =  anc*  t^le  integral  in  the  denom¬ 
inator  of  (7)  is  finite.  Thus  the  updated  pdf  Fk  v  after  the 
iF  consensus  loop  is  the  weighted  geometric  average  of  the 
pdfs  of  the  inclusive  neighbors  Fkv_kfil  £  from  the 
{v  —  l)th  consensus  loop,  at  the  AJh  time  instant.  The  LogOP 
solution  is  typically  unimodal  and  less  dispersed,  indicating  a 
consensual  estimate  jointly  preferred  by  the  network  [25].  The 
LogOP  solution  is  invariant  under  under  rescaling  of  individual 
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degrees  of  belief,  hence  it  preserves  an  important  credo  of  uni- 
Bayesian  decision  theory;  i.e.,  the  optimal  decision  should  not 
depend  upon  the  choice  of  scale  for  the  utility  function  or 
prior  probability  distribution.  The  most  compelling  reason  for 
using  LogOP  is  that  it  is  externally  Bayesian;  i.e.,  finding  the 
consensus  distribution  commutes  with  the  process  of  revising 
distributions  using  a  commonly  agreed  likelihood  distribution. 
Next,  we  present  consensus  theorems  using  the  LogOP. 

Assumption  5.  The  local  estimated  pdf  at  the  start  of  the  con¬ 
sensus  stage  is  positive  everywhere,  i.e.,  Fk  0  =  pk0(xk)  > 
0,  Vxk  £  X,\/j  £  □ 

Assumption  5  is  introduced  to  avoid  regions  with  zero  prob¬ 
ability,  since  they  would  constitute  vetoes  and  unduly  great 
emphasis  would  get  placed  on  them.  Moreover,  the  LogOP 
guarantees  that  Pk  v  will  remain  positive  for  all  subsequent 
consensus  loop. 

Definition  1.  (Hi  vector  for  LogOP)  For  the  purpose  of 
analysis,  let  us  choose  x^o  £  X  such  that  pk  „(xko)  >  0,  Vj  € 

{1 , . . .  €  N.  Let  us  define  Hk  v  :=  In 

Under  Assumption  5,  Hk  v  is  a  well-defined  function,  but  need 
not  be  a  C  \  function.  Then,  by  simple  algebraic  manipulation 
of  (7),  we  get; 

nk,u=  E  €  {L,---,rn},v  £  N,  (8) 

ieJi 

which  is  similar  to  the  LinOP  (6).  Let  Uk  v  '■= 

\  T 

H\  v, . . . ,  H™v  j  be  an  array  of  the  estimates  of  all  the 
agents  during  the  ;/lh  consensus  loop  at  the  /,:lh  time  instant, 
then  the  equation  (8)  can  be  expressed  concisely  as: 


PkAxl  o) 


•F k  ““  PkiXk )  — 


n?=i(pj,0(*k))” 

IxUZi  (pi,o(xk))m  Mx) 


(10) 


at  a  rate  faster  or  equal  to  Am_i (Pk  Pk)  =  crm_i(Pfc). 
Furthermore,  their  induced  measures  globally  exponentially 

TV. 

converge  in  total  variation,  i.e.,  lim,,-^  pTj  -A-  /ijr* ,  V)  £ 

to}. 


Proof:  See  Appendix  A.  ■ 

The  KL  divergence  is  a  measure  of  the  information  lost 
when  the  consensual  pdf  is  used  to  approximate  the  locally 
estimated  posterior  pdfs.  We  now  show  that  the  consensual 
pdf  Tk  obtained  using  Theorem  2,  which  is  the  weighted 
geometric  average  of  the  locally  estimated  posterior  pdfs 
Tk  o .  Vj  £  {1, . . . ,  m},  minimizes  the  information  lost  during 
the  consensus  stage  because  it  minimizes  the  sum  of  KL 
divergences  with  those  pdfs. 


Theorem  3.  The  consensual  pdf  Tk  given  by  (10)  globally 
minimizes  the  sum  of  Kullback-Leibler  (KL)  divergences  with 
the  locally  estimated  posterior  pdfs  at  the  start  of  the  consen¬ 
sus  stage  (FJk0,\/j  £  {1, . . . ,  m},  i.e., 

m 

J-*  =  arg  min  ^  dkl  (p\  \H,o)  -  (u) 

pGJt ?i(X)  “ 

t=i 

where  Jfi(X)  is  the  set  of  all  pdfs  over  the  state  space  X 
satisfying  Assumption  5. 


Proof:  The  sum  of  the  KL  divergences  of  a  pdf  p  £ 
Jfi(X)  with  the  locally  estimated  posterior  pdfs  is  given  by; 


Uk,v  =  Pk,v- \Uk.v-i,  VP  G  N,  (9) 

where  Pk,v- 1  is  a  matrix  with  entries  ak  [/_1.  □ 

Thus  we  are  able  to  use  the  highly  nonlinear  LogOP 
for  combining  the  pdf  estimates,  but  we  have  reduced  the 
complexity  of  the  problem  to  that  of  consensus  using  the 
LinOP.  Next,  we  discuss  the  algorithm  for  achieving  global 
exponential  convergence  on  balanced  graphs  using  the  LogOP. 

Assumption  6.  The  communication  network  topology  of 
the  multi-agent  system  Qk  is  strongly  connected  (SC)  and 
balanced.  The  weighting  factors  akiy_1,Wj,£  £  {1, 
and  the  matrix  Pk,v- 1  have  the  following  properties;  (i)  the 
weighting  factors  are  the  same  for  all  consensus  loops  within 
each  time  instant,  (ii)  the  matrix  Pj,  conforms  with  the  graph 
Qk,  (iii)  the  matrix  Pk  is  row  stochastic,  and  (iv)  the  weighting 
factors  a:k  are  such  that  the  digraph  Qk  is  balanced.  □ 

Theorem  2.  (Consensus  using  the  LogOP  on  SC  Balanced 
Digraphs)  Under  Assumption  5  and  6,  using  the  LogOP  (7), 
each  (Fk  globally  exponentially  converges  pointwise  to  the 
Pdf  Pk  given  by: 


m 

(p\\H,o)  = 

i= 1 

m  r. 

V  /  (p(xk')Mp(xk))-p(xk)ln(plkfi(xk)))  dp(x).  (12) 

r=lJx 

Taking  the  logarithm,  in  the  KL  divergence  formula,  to  the 
base  c  :=  (^$XY\Z=\  {Pk,o(xk)^j  ™  dp(x)^j  and  then  differen¬ 
tiating  YZZ i  Dkl  {p\\Plk  0)  with  respect  to  p  using  Leibniz 
integral  rule  [36,  pp.  372]  gives; 

m  r. 

E  /  (loSc (p(*fc))  +  1  -  logc(Pfc)0(*fc)))  dp(x)  =  0  , 

i=l Jx 

which  is  minimized  by  Pk.  ■ 

Note  that  if  a  central  agent  receives  all  the  locally  estimated 
posterior  pdfs  (Pk  Q;  Vj  £  {l,...,m})  and  is  tasked  to  find 
the  best  estimate  in  the  information  theoretic  sense,  then  it 
would  also  yield  the  same  consensual  pdf  Fk  given  by  (10). 
Hence  we  claim  to  have  achieved  distributed  estimation  using 
this  algorithm. 
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B.  Communicating  Probability  Distributions 

The  consensus  algorithms  need  the  estimated  pdfs  to  be 
communicated  to  other  agents  in  the  network.  The  first  ap¬ 
proach  involves  approximating  the  pdf  by  a  weighted  sum  of 
Gaussians  and  then  transmitting  this  approximate  distribution. 
Let  A f(xk  —  rri, .  II, )  denote  the  Gaussian  density  function, 
where  the  mean  is  the  nx-vector  ru,  and  the  covariance  is 
the  positive -definite  symmetric  matrix  P*.  The  Gaussian  sum 
approximations  lemma  of  [38,  pp.  213]  states  that  any  pdf 
F  —  p(xk)  can  be  approximated  as  closely  as  desired  in 
the  2z?i(lR”*)  space  by  a  pdf  of  the  form  F  =  p(xk)  = 
Y^i=iai N{xk  —  ),  for  some  integer  ng  and  positive 

scalars  a,;  with  X]”=i  ai  =  1-  For  an  acceptable  communica¬ 
tion  error  ecomm  >  0,  there  exists  ng,  ai,  rrii  and  Bi  such  that 
\\?  -  F\\Cl  <  ecomm.  Let  Fi  be  the  LinOP  solution  after 
combining  local  pdfs  corrupted  by  communication  error,  i.e., 
F:fv  :=  T  where  T(-)  is  either  LinOP  or 

LogOP.  Then  we  get  \\Pk  v  -  Pk  v <  ve  comm?  Mv  e  N, 
where  Fi  is  the  true  solution  obtained  from  uncormpted 
local  pdfs.  Several  techniques  for  estimating  the  parameters 
are  discussed  in  the  Gaussian  mixture  model  literature,  like 
maximum  likelihood  (ML)  and  maximum  a  posteriori  (MAP) 
parameter  estimation  [39]-  [41],  Hence,  in  order  to  commu¬ 
nicate  the  pdf  F,  the  agent  needs  to  transmit  \ngnx  (nx  +  3) 
real  numbers. 

If  particle  filters  are  used  to  evaluate  the  Bayesian  filter  and 
combine  the  pdfs  [22],  then  the  resampled  particles  represent 
the  agent’s  estimated  pdf  of  the  target.  Hence  communicating 
pdfs  is  equivalent  to  transmitting  these  resampled  particles. 
The  information  theoretic  approach  for  communicating  pdfs 
is  discussed  in  [42],  Now  that  we  have  established  that 
communication  of  pdfs  is  possible,  let  us  discuss  the  complete 
BCF  algorithm. 


Algorithm  1  BCF-LogOP  on  SC  Balanced  Digraphs 

1 

(one  cycle  of  jth  agent  during  fc,h  time  instant) 

2 

Given  the  pdf  from  previous  time  ste 

P 

H- 1  =Pfc-i(*k-i) 

3 

Set  nioop,  the  weights  aijf  } 

Theorems  2,  4 

4 

while  tracking  do 

5 

Compute  the  prior  pdf  pJk(xk) 

Bayesian 

6 

Compute  the  posterior  pdf  Fj 

Filtering  Stage 

7 

for  v  =  1  to  nioop 

8 

if  v  =  1  then  Set  F3k  0  =  F3k 

end  if 

LogOP-based 

9 

Obtain  the  communicated 

>  Consensus 

Pdfs  Pfc,P_i,V£  €  Jl 

Stage 

10 

Compute  the  new  pdf  Fk  v 

end  for 

11 

Set  Fl  =  Fl  end  while 

O'  O' )  '  Hoop 

IV.  Main  Algorithms:  Bayesian  Consensus 
Filtering 

In  this  section,  we  finally  present  the  BCF  algorithm  il¬ 
lustrated  by  Algorithm  1.  The  BCF  is  performed  in  two 


steps:  (i)  each  agent  locally  estimates  the  pdf  of  the  target’s 
states  using  a  Bayesian  filter  with/without  measurements  from 
neighboring  agents,  as  discussed  in  Section  II-A,  and  (ii) 
during  the  consensus  stage,  each  agent  recursively  transmits 
its  pdf  estimate  of  the  target’s  states  to  other  agents,  receives 
estimates  of  its  neighboring  agents,  and  combines  them  using 
the  LogOP  as  discussed  in  Section  III-A.  In  this  section,  we 
compute  the  number  of  consensus  loops  (nioop)  needed  to 
reach  a  satisfactory  consensus  estimate  across  the  network  and 
discuss  the  convergence  of  this  algorithm. 

Definition  2.  ( Disagreement  vector  9klf)  Let  us  define 

ekM  :=  0£„)T,  Where  9{^  :=  \\F^  -  F£\\Cl. 

Since  the  C\  distance  between  pdfs  is  upper  bounded  by  2, 
the  li  norm  of  the  disagreement  vector  (||0fc,i'||r2)  is  upper 
bounded  by  2  i/to.  □ 

This  conservative  bound  is  used  to  obtain  the  minimum 
number  of  consensus  loops  for  achieving  e-consensus  across 
the  network,  while  tracking  a  moving  target.  Let  us  now 
quantify  the  divergence  of  the  local  pdfs  during  the  Bayesian 
filtering  stage. 

Definition  3.  ( Error  propagation  dynamics  T(-))  Let  us  as¬ 
sume  that  the  dynamics  of  the  £ 2  norm  of  the  disagreement 
vector  during  the  Bayesian  filtering  stage  can  be  obtained  from 
the  target  dynamics  and  measurement  models  (1)  and  (2). 
The  error  propagation  dynamics  T(-)  estimates  the  maximum 
divergence  of  the  local  pdfs  during  the  Bayesian  filtering  stage, 
i-e.,  ||0fc,o|k  <  r  (||0fc-i,„loop||<?2),  where  ||6>fe-i,moop |U2  is  the 
disagreement  vector  with  respect  to  Fk_\  at  the  end  of  the 
consensus  stage  during  the  (k  —  l)th  time  instant;  and  ||0fc,o||r2 
is  the  disagreement  vector  with  respect  to  Fk  after  the  update 
stage  during  the  fc*  time  instant.  □ 

Next  we  obtain  the  minimum  number  of  consensus  loops 
for  achieving  ^-consensus  across  the  network  and  also  derive 
conditions  on  the  communication  network  topology  for  a  given 
number  of  consensus  loops. 

Theorem  4.  (BCF-LogOP  on  SC  Balanced  Digraphs)  Under 
Assumptions  5,  6,  and  in  the  absence  of  communication 
inaccuracies,  each  agent  tracks  the  target  using  the  BCF  al¬ 
gorithm.  For  some  acceptable  consensus  error  ^consensus  >  0 
and  -ik  =  min  (r  (||0fc_i,„loop  \\e2)  ,2 s/m): 

(i)  if  the  number  of  consensus  loops  is  at  least  ?ii00p  > 

f°r  a  given  Pk;  or 

(ii)  if  the  communication  network  topology  (Pk)  during  the  kth 
time  instant  is  such  that  am-i(Pk)  <  ^ ^otaemu. ^  loop  yor  a 
given  nioop/ 

then  the  £ 2  norm  of  the  disagreement  vector  at  the  end  of  the 
consensus  stage  is  less  than  £consensVlS,  i.e.,  ||0fc,moop  \\l2  < 

^consensus- 

Proof:  It  follows  from  Theorem  2  that,  if  9k  o  is 
the  initial  disagreement  vector  at  the  start  of  the  consen¬ 
sus  stage,  then  ||0fc>nioop \\e2  <  (frm_i(Pfc))"loop  ||0fc,o|k  < 
(trm_i(Pfc))riloop  7fc.  Thus,  we  get  the  conditions  on  nioop 
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or  <7m_i(Pfe)  from  the  inequality  (crm_i(Pfc))nioop  7fc  < 

^consensus*  ® 

Note  that  for  the  particular  case  where  nioop  =  1,  we  need 
<rm_i(Pfc)  <  ' coiiscnaua  for  g-convergence  across  the  network. 


A.  Hierarchical  Bayesian  Consensus  Filtering 

In  this  section,  we  modify  the  original  problem  statement 
such  that  only  Toi  out  of  to  agents  are  able  to  observe 
the  target  at  the  fcth  time  instant.  In  this  scenario,  the  other 
m2  (=  to  —  toi)  agents  are  not  able  to  observe  the  target. 
Without  loss  of  generality,  we  assume  that  the  first  toi 
agents,  i.e.,  {1,  2, . . . ,  Toi},  are  tracking  the  target.  During 
the  Bayesian  filtering  stage,  each  tracking  agent  (i.e.,  agent 
tracking  the  target)  estimates  the  posterior  pdf  of  the  target’s 

states  at  the  fc*  time  instant  {PJk  =  jPk(xk\z^,k^1'"'mi^ ) ,  Vj  £ 
{ 1 , . . . ,  Toi })  using  the  estimated  prior  pdf  of  the  target’s 

states  and  the  new  measurement  array  zk  := 

jz^.,  W  £  S3k  n  {1, . . .  TOi}|  obtained  from  the  neighboring 
tracking  agents.  Each  non-tracking  agent  (i.e.,  agent  not  track¬ 
ing  the  target)  only  propagates  its  prior  pdf  during  this  stage 
to  obtain  pjk(xk),\/j  £  {toi  +  1, . . .  ,to}. 

The  objective  of  hierarchical  consensus  algorithm  is  to  guar¬ 
antee  pointwise  convergence  of  each  7Fk  £  {1, . . .  ,to} 
to  a  pdf  JFk  and  only  the  local  estimates  of  the  agents  tracking 
the  target  contribute  to  the  consensual  pdf.  This  is  achieved  by 
each  tracking  agent  recursively  transmitting  its  estimate  of  the 
target’s  states  to  other  agents,  only  receiving  estimates  from  its 
neighboring  tracking  agents  and  updating  its  estimate  of  the 
target.  On  the  other  hand,  each  non-tracking  agent  recursively 
transmits  its  estimate  of  the  target’s  states  to  other  agents, 
receives  estimates  from  all  its  neighboring  agents  and  updates 
its  estimate  of  the  target.  Let  T>k  represent  the  communication 
network  topology  of  only  the  tracking  agents. 

Assumption  7.  The  communication  network  topologies  Qk 
and  T>k  are  SC  and  T)k  is  balanced.  The  weighting  factors 
«i£„_i>Vj, £  £  {l,...,m}  and  the  matrix  Pk,v-\  have  the 
following  properties:  (i)  the  weighting  factors  are  the  same 
for  all  consensus  loops  within  each  time  instants;  (ii)  if  j  £ 
{1, . . . ,  toi},  then  ak  >  0  if  and  only  if  i  £  J^C |{1, . . . ,  toi}, 
else  ak  =  0;  (iii)  if  j  £  {toi  +  1, . . . ,  to},  then  aJk  >  0  if 
and  only  if  f  £  J7j?,  else  aJk  =  0;  (iv)  the  matrix  Pk  is  row 
stochastic,  and  (v)  the  weighting  factors  aJk  are  such  that  the 
digraph  T>k  is  balanced.  □ 


Theorem  5.  (Hierarchical  Consensus  using  the  LogOP  on 
SC  Balanced  Digraphs)  Under  Assumptions  5  and  7,  using 
the  LogOP  (7),  each  Tk  v  globally  exponentially  converges 
pointwise  to  the  pdf  Tk  given  by: 


H  =  Pk(xk )  = 


uz\  (pi,o(^)j 


fx  nr=i  (pm (**))  "ll  Mx) 


(13) 


at  a  rate  faster  or  equal  to  y  Ami_i(P^}P/ci)  =  <rmi_i(Pfci), 
where  the  matrix  Pk\  conforms  with  the  directed  graph  T>k. 


Fig.  1.  The  SSN  locations  are  shown  along  with  their  static  SC  balanced 
communication  network  topology.  The  orbit  of  the  Iridium-33  debris  is  shown 
in  red,  where  marks  its  actual  position  during  particular  time  instants. 


Only  the  initial  estimates  of  the  tracking  agents  contribute  to 
the  consensual  pdf  Tk.  Furthermore,  their  induced  measures 

TV. 

converge  in  total  variation,  i.e.,  lim,,-^  p^j  -A-  p,jr* ,  Vj  £ 

{1, . . . ,  to}. 

Proof:  The  proof  follows  from  Theorem  2.  In  fact,  the 
inessential  states  die  out  geometrically  fast  [43,  Theorem  4.3, 

pp.  120].  ■ 

A  simulation  example  of  Hierarchical  BCF-LogOP  algo¬ 
rithm  for  tracking  orbital  debris  in  space  is  discussed  in  the 
next  section. 


V.  Numerical  Example 

Currently,  there  are  over  ten  thousand  objects  in  Earth 
orbit,  of  size  0.5  cm  or  greater,  and  almost  95%  of  them  are 
nonfunctional  space  debris.  These  debris  pose  a  significant 
threat  to  functional  spacecraft  and  satellites  in  orbit.  The  US 
has  established  the  Space  Surveillance  Network  (SSN)  for 
ground  based  observations  of  the  orbital  debris  using  radars 
and  optical  telescopes  [44]  (See  Fig.  1).  In  February  2009,  the 
Iridium-33  satellite  collided  with  the  Kosmos-2251  satellite 
and  a  large  number  of  debris  fragments  were  created.  In  this 
section,  we  use  the  Hierarchical  BCF-LogOP  Algorithm  to 
track  one  of  the  Iridium-33  debris  created  in  this  collision. 

The  actual  two-line  element  set  (TLE)  of  the  Iridium- 
33  debris  was  accessed  from  North  American  Aerospace 
Defense  Command  (NORAD)  on  4th  Dec  2013.  The  nonlinear 
Simplified  General  Perturbations  (SGP4)  model,  which  uses  an 
extensive  gravitational  model  and  accounts  for  the  drag  effect 
on  mean  motion  [45],  [46],  is  used  as  the  target  dynamics 
model.  If  the  debris  is  visible  above  the  sensor’s  horizon,  then 
it  is  assumed  to  create  a  single  measurement  during  each  time 
step  of  one  minute.  The  heterogeneous  measurement  model  of 
the  jth  sensor  is  given  by: 

zJk  =  xk  +  w\ where  wJk  =  A f  (0,  (1000  +  50j)  x  I) , 

where  xk  £  R3  is  the  actual  location  of  the  debris.  Since  it 
is  not  possible  to  implement  the  SGP4  target  dynamics  on 
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Fig.  2.  (a)  Number  of  SSN  sensors  observing  debris,  (b)  Trajectories  of 

particles  for  stand-alone  Bayesian  filters  for  10th  SSN  sensor.  Trajectories  of 
particles  of  all  sensors  for  (c)  Hierarchical  BCF-LinOP  and  (d)  Hierarchical 
BCF— LogOP. 

distributed  estimation  algorithms  discussed  in  the  literature 
[1]-  [15],  we  compare  the  performance  of  our  Hierarchical 
BCF-LogOP  algorithm  against  the  Hierarchical  BCF-LinOP 
algorithm,  where  the  LinOP  is  used  during  the  consensus 
stage. 

In  this  simulation  example,  we  simplify  the  debris  tracking 
problem  by  assuming  only  the  mean  motion  (n)  of  the  debris 
is  unknown,  which  needs  to  be  estimated  within  100  minutes. 
Hence,  each  sensor  knows  the  other  TLE  parameters  of  the 
debris  and  an  uniform  prior  distribution  (Pq)  is  assumed.  Note 
that  at  any  time  instant,  only  a  few  of  the  SSN  sensors  can 
observe  the  debris,  as  shown  in  Fig  2(a).  The  estimates  of  the 
stand-alone  Bayesian  filter  for  the  10th  sensor  do  not  converge 
due  to  large  measurement  errors,  in  spite  of  observing  the 
debris  for  some  time. 

Particle  filters  with  resampling  are  used  to  evaluate  the 
Bayesian  filters  and  communicate  pdfs  in  the  Hierarchical 
BCF  algorithms.  100  particles  are  used  by  each  sensor  and 
10  consensus  loops  are  executed  during  each  time  step  of  one 
minute.  As  expected,  all  the  sensors  converge  on  the  correct 
value  of  n  of  14.6  revs  per  day.  The  Hierarchical  BCF-LinOP 
estimates  are  multimodal  for  the  first  90  minutes.  On  the 
other  hand,  the  Hierarchical  BCF-LogOP  estimates  converges 
to  the  correct  value  within  the  first  10  minutes  because  the 
LogOP  algorithm  efficiently  communicates  the  best  consensual 
estimate  to  other  sensors  during  each  time  step  and  achieves 
consensus  across  the  network. 

VI.  Conclusion 

In  this  paper,  we  extended  the  scope  of  distributed  esti¬ 
mation  algorithms  in  a  Bayesian  filtering  framework  in  order 
to  simultaneously  track  targets,  with  general  nonlinear  time- 
varying  target  dynamic  models,  using  a  strongly  connected 


network  of  heterogeneous  agents,  with  general  nonlinear  time- 
varying  measurement  models.  The  LogOP  algorithm  on  SC 
balanced  digraph  converges  globally  exponentially,  and  the 
consensual  pdf  minimizes  the  information  lost  during  the 
consensus  stage  because  it  minimizes  the  sum  of  KL  di¬ 
vergences  to  each  locally  estimated  probability  distribution. 
We  introduced  the  BCF  algorithm,  where  the  local  estimated 
posterior  pdfs  of  the  target’s  states  are  first  updated  using 
the  Bayesian  filter  and  then  recursively  combined  during  the 
consensus  stage  using  LogOP,  so  that  the  agents  can  track  a 
moving  target  and  also  maintain  consensus  across  the  network. 
Conditions  for  exponential  convergence  of  the  BCF  algorithm 
and  constraints  on  the  communication  network  topology  have 
been  studied.  The  Hierarchical  BCF  algorithm,  where  some  of 
the  agents  do  not  observe  the  target,  has  also  been  investigated. 
Simulation  results  demonstrate  the  effectiveness  of  the  BCF 
algorithms  for  nonlinear  distributed  estimation  problems. 


Appendix  A 
Proof  of  Theorem  2 

Under  Assumption  6,  Pk  is  a  nonnegative,  doubly  stochastic 
and  irreducible  matrix.  Hence  Perron-Frobenius  theorem  (cf. 
[43,  pp.  3])  states  that  lim^o 0  Pk  =  ^llr  and  each  H3kv 
converges  pointwise  to  ± 1 TUkfi  =  ±  Y^Li^-1,0-  We 

have  \/xk  €  X: 

hm  (ln^  ^-lnpi^o))  =  lnp*k(xk)-hip*k(xk0). 

As  3xk0  £  X  such  that  lim^oc  p{^(,xk0)  =  pl(xk0),  we 
get  lim^o cp{^(xk)  =  pl{xk),Mxk  £  P.  Thus  each  T3kv 
converges  pointwise  to  the  consensual  pdf  P;*  given  by  (10). 
By  Lemma  1,  the  measure  induced  by  Tk  v  on  3P  converges 

in  total  variation  to  the  measure  induced  by  T*k  on  ,P,  i.e., 
,.  T.V. 

lim^oo  M  to  - >  fP;  ■ 


If  Vtr  = 


4=1,  Us 

,  /TO.  '  a 


are  the  orthonormal  eigenvectors  of 


Pk  Pk,  then  by  spectral  decomposition  [47]  we  get  that  the 
synchronizes  to  -^=1  (or  Uk)  is  equal 


y/rn 

0(m— l)xl.  If  $  = 


rate  at  which  Uk 

to  the  rate  at  which  V,T/4, „  —t  u'  n  'vk,v 

(y^UktV)TV^Uk^  is  a  candidate  Lyapunov  function,  then 
<  (Ama x(Us^PjPfcUs))  l-  Hence  each  Hkv  glob¬ 

ally  exponentially  converges  pointwise  to  1~Lk  with  a  rate  faster 

or  equal  to  \m_i{P£ Pk)  =  am-i(Pk). 

Next,  we  need  to  find  the  rate  of  convergence  of  P,' 
to  Tk.  Let  us  define  the  continuous  function  al  (xk)  such 


that  aJk:Jxk )  = 


Pk,  v(Xk)  P*kixkO ) 

Pl(Xk)  pJk  ^(Xk0) 


if  Pk,u(Xk)p*k{xk0)  > 


Pk(xk)pi,v(xkO )  and  a3k  u(xk)  = 

wise.  Then  we  get: 


Pk(xh)  pj.JXk  o) 

PJkJxk)  PfcOfco) 


other- 


a) <  («fc,o(*fc)) 


i  {Pk)Y 


(14) 


Using  the  mean  value  theorem,  (14)  can  be  simplified  to: 


aiAxk)  - 1  <  ( <jm-i{Pk)Y  (o£,0 {xk) 


- 1 


(15) 
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Irrespective  of  the  orientation  of  aj,  v{xk)  and  ctJk  o(xk)< 
(15)  can  be  written  as  (16)  by  multiplying  with  ;  1 — -  or 

is  &  / 

'{xiy  and  then  with  p*k(xk). 


PUxkO ) 
Pk,AXk0) 


Pk,Axk)  ~  P*k(xk) 


<  ((Jm-l(Pk))1' 


Pk(.xkO ) 
p]k0(xk0) 


pI,o  (xk)~P*k(xk) 


(16) 


If  we  choose  Xko  €  X  such  that  p^.0(&feo)  =  Pk(xk o),  then 
we  can  simplify  (16)  to: 


Pk,v(xk)  ~Pk(xk) 


<  ( am-i{Pk)Y 


Pk,0(Xk)  -  Pk(Xk ) 


Thus  each  Fi  u  =  pi  „(xk)  globally  exponentially  converges 
to  Fk  =  p*k{xk)  with  a  rate  faster  or  equal  to  <jm_i(Pfc).  ■ 
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